library(dplyr)
library(MASS)
library(candisc)
library(klaR)
library(FSA)
library(ROCR)
library("PerformanceAnalytics")
library("FactoMineR")
library("factoextra")
Возьмем набор данных Wine quality dataset.
Этот набор данных включает данные о физико-химическом составе красных вин.
df_raw <- read.csv("winequality-white.csv", sep = ";")
head(df_raw)
Признаки:
fixed.acidity - фиксированная кислотность.volatile.acidity- летучая кислотность.citric.acid - лимонная кислота.residual.sugar - остаточный сахар.chlorides - хлориды.free.sulfur.dioxide - свободный диоксид серы.total.sulfur.dioxide - общий диоксид серы.density - плотность.pH - pH.sulphates - сульфаты.alcohol - алкоголь.quality - качество.Количество индивидов:
length(df_raw[, 1])
## [1] 4898
Посмотрим на количество классов.
count(df_raw, quality)
И на их соотношение:
perc_class <- count(df_raw, quality)
perc_class$n <- count(df_raw, quality)$n / sum(count(df_raw, quality)$n)
perc_class
Обозначим классы qualuity:
1 (“ниже среднего”)2 (“выше среднего”)3 (“отличные”)df <- df_raw %>% filter(quality > 4, quality < 8)
df_low <- df_raw %>% filter(quality < 5)
df_low$quality <- 5
df_high <- df_raw %>% filter(quality > 7)
df_high$quality <- 7
df <- rbind(df, df_low, df_high)
df$quality <- cut(df$quality, 3, labels=c('1', '2', '3'))
Посмотрим, есть ли признаки, которые стоит логарифмировать.
chart.Correlation(subset(df, select=-c(quality)), histogram=TRUE)
Прологарифмируем residual.sugar, volatile.acidity, citric.acid, chlorides, free.sulfur.dioxide, sulphates, density.
df_l <- df
df_l <- transform(df_l, residual.sugar=log(residual.sugar))
df_l <- transform(df_l, chlorides=log(chlorides))
df_l <- transform(df_l, citric.acid=log(citric.acid))
df_l <- transform(df_l, volatile.acidity=log(volatile.acidity))
df_l <- transform(df_l, free.sulfur.dioxide=log(free.sulfur.dioxide))
df_l <- transform(df_l, sulphates=log(sulphates))
df_l <- transform(df_l, density=log(density))
chart.Correlation(subset(df_l, select=-c(quality)), histogram=TRUE)
Есть 19 индивидов, у которых citric.acid был равен нулю. Удалим эти индивиды.
df_l <- df_l %>% filter(!is.infinite(citric.acid))
Разобьем выборку на \(train\) и \(test\).
smp_size <- floor(0.75 * nrow(df))
set.seed(123)
train_ind <- sample(seq_len(nrow(df_l)), size = smp_size)
train <- df_l[train_ind, ]
test <- df_l[-train_ind, ]
rownames(train) = 1:length(train[, 1])
rownames(test) = 1:length(test[, 3])
Посмотрим на баланс классов в \(train\) и \(test\).
count(train, quality)
count(test, quality)
И на их процентное соотношение:
train:perc_class <- count(train, quality)
perc_class$n <- count(train, quality)$n / sum(count(train, quality)$n)
perc_class
test:perc_class <- count(test, quality)
perc_class$n <- count(test, quality)$n / sum(count(test, quality)$n)
perc_class
Проверим гипотезу о равенстве мат. ожиданий среди групп quality.
train.manova <- manova(cbind(fixed.acidity, volatile.acidity, citric.acid, residual.sugar, chlorides, free.sulfur.dioxide, total.sulfur.dioxide, density, pH, sulphates, alcohol) ~ quality, data = train)
summary(train.manova, 'Wilks')
## Df Wilks approx F num Df den Df Pr(>F)
## quality 2 0.67848 71.217 22 7320 < 2.2e-16 ***
## Residuals 3670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(train.manova, 'Roy')
## Df Roy approx F num Df den Df Pr(>F)
## quality 2 0.44341 147.58 11 3661 < 2.2e-16 ***
## Residuals 3670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Однако перед тем, как интерпретировать результаты критериев, посмотрим еще раз на распределение прологарифмированных данных. Попробуем понять хотя бы приблизительно, равны ли ковариационные матрицы среди групп quality.
ggplot(df_l) +
geom_point(aes(x = pH, y = alcohol, group = quality, color = quality)) +
stat_ellipse(aes(x = pH, y = alcohol, group = quality, color = quality)) +
scale_color_manual(values = c("red", "purple", "black"))
ggplot(df_l) +
geom_point(aes(x = fixed.acidity, y = volatile.acidity, group = quality, color = quality)) +
stat_ellipse(aes(x = fixed.acidity, y = volatile.acidity, group = quality, color = quality)) +
scale_color_manual(values = c("red", "purple", "black"))
По формам эллипсов можно уже заметить, что ковариационные матрицы различны, откуда следует тот факт, что на критерии Wilk’s Lambda или Roy’s greatest root нельзя полагаться.
Более того, из определения модели \(LDA\) плотность в точке \(x\): \(p_i(x) = p(x|\xi = A_i) = \frac{1}{2\pi^{p/2}|\Sigma|^{1/2}}exp(-\frac{1}{2}(x-\mu_i)^\mathrm{T}\Sigma^{-1}(x-\mu_i))\) и классифицирующие функции \(f_i(x) = \pi_ip(x|\xi = A_i)\) предполагают равенство ковариационных матриц между классами. Хотя мы и проверяем “на глаз”, скорее всего нам больше подойдет \(QDA\), нежели \(LDA\).
Обучим модель на train и на нём же протестируем.
full.lda <- lda(x=train[, 1:11], grouping = train[, 12])
full.ldap <- predict(full.lda, train[, 1:11])
full.ldapc <- full.ldap$class
ct <- table(full.ldapc, train[, 12])
ct
##
## full.ldapc 1 2 3
## 1 700 332 30
## 2 491 1141 446
## 3 21 200 312
Точность модели по классам:
diag(prop.table(ct, 2))
## 1 2 3
## 0.5775578 0.6820084 0.3959391
Попробуем задать априорные вероятности \(\pi_i\) равными, а не на основе соотношения количества элементов в выборке.
full.lda <- lda(x=train[, 1:11], grouping = train[, 12], prior = c(1, 1, 1)/3)
full.ldap <- predict(full.lda, train[, 1:11])
full.ldapc <- full.ldap$class
ct <- table(full.ldapc, train[, 12])
ct
##
## full.ldapc 1 2 3
## 1 857 502 85
## 2 257 599 154
## 3 98 572 549
Точность (accuracy) модели по классам:
diag(prop.table(ct, 2))
## 1 2 3
## 0.7070957 0.3580395 0.6967005
Видим что точность для \(1\)го и \(3\)го класса сильно возрасла, в то время как для \(2\)го класса упала.
Посмотрим на результаты при обучении на train и проверке на test.
Обучим модель:
train.lda <- lda(train[, 1:11], train[, 12], prior = c(1, 1, 1)/3)
train.ldap <- predict(train.lda, test[, 1:11])
train.ldapc <- train.ldap$class
ct <- table(train.ldapc, test[,12])
ct
##
## train.ldapc 1 2 3
## 1 291 160 27
## 2 89 168 56
## 3 37 189 189
Точность модели по классам:
diag(prop.table(ct, 2))
## 1 2 3
## 0.6978417 0.3249516 0.6948529
Построим модель используя кросс-валидацию:
cv.lda <- lda(df_l[, 1:11], df_l[, 12], CV = TRUE, prior = c(1, 1, 1)/3)
cv.ldac <- cv.lda$class
ct <- table(cv.ldac, df_l[,12])
ct
##
## cv.ldac 1 2 3
## 1 1170 670 116
## 2 336 757 217
## 3 123 763 727
diag(prop.table(ct, 2))
## 1 2 3
## 0.7182320 0.3456621 0.6858491
result.candisc <- candisc(train.manova)
ggplot(result.candisc$scores) +
geom_point(aes(x = Can1, y = Can2, color = quality)) +
stat_ellipse(aes(x = Can1, y = Can2, color = quality)) +
geom_segment(as.data.frame(result.candisc$structure),
mapping = aes(x = 0, y = 0, xend = 6.5 * Can1, yend = 6.5 * Can2),
arrow = arrow(angle = 10, type = "closed"), size = 0.5)
Изображение разделяющих плоскостей в разных проекциях
plot(partimat(quality ~ ., train, method="lda"), xlim = c(-10, 10),
ylim = c(-10, 10));